#ifndef EQUATIONSOLVER
#define EQUATIONSOLVER

#include "Function.hpp"
#include <iostream>
#include <cmath>
#include <stdexcept>
class EquationSolver{
protected:
    const Function & F;
public:
    EquationSolver(const Function& F) : F(F) {}
    virtual double solve() = 0;
};

class Bisection_Method : public EquationSolver {
private:
    double a, b;
    double eps, delta;
    int Maxiter;
public:
    Bisection_Method(const Function &F, double a, double b, 
        double eps = 1e-7, double delta = 1e-6, int Maxiter = 50) :
        EquationSolver(F), a(a), b(b), eps(eps), delta(delta), Maxiter(Maxiter) {
		// check whether the input valid
        if (F(a) * F(b) > 0) {
            throw std::invalid_argument("The function values at the interval endpoints must have opposite signs.");
        }
		}
    
    virtual double solve() {
    	
    	double c = (b - a)/2;
        for(int i = 0; i < Maxiter; i++){
        	if(std::abs(F(a + c)) < eps || c < delta){
        		return a + c;
			}
		if(F(a)*F(a +c) > 0){
			a = a + c;
		}
		c  = c/2;
		}
    }
};

class Newton_Method : public EquationSolver {
private:
    double x0;
    double eps;
    int Maxiter;
public:
    Newton_Method(const Function &F, double x0, 
        double eps = 1e-7, double Maxiter = 10) :
        EquationSolver(F), x0(x0), Maxiter(Maxiter), eps(eps) {}
    
    virtual double solve() {
        double x1 = x0;
	    for (int i = 0; i < Maxiter; i++) {
	    	if (std::abs(F(x1)) < eps){
	    		return x1; 
			} 
	        double fx = F(x1);
	        double dfx = F.derivative(x1);
	        if (dfx == 0) {
	            std::cerr << "Derivative is zero, no solution found." << std::endl;
	            return x1;
	        }
	        x1 = x1 - fx / dfx;
	    }
	    return x1;
    }
};

class Secant_Method : public EquationSolver {
private:
    double x0, x1;
    double eps, delta;
    int Maxiter;
public:
    Secant_Method(const Function &F, double x0, double x1, 
        double eps = 1e-7, double delta = 1e-6, double Maxiter = 15) :
        EquationSolver(F), x0(x0), x1(x1), Maxiter(Maxiter), eps(eps), delta(delta) {}
    
    virtual double solve() {
	    for (int i = 0; i < Maxiter; i++) {
	    	double fx = F(x1);
	    	if (std::abs(fx) < eps || std::abs(x1 - x0) < delta){
	    		return x1; 
			} 
	    	double s = (x1 - x0)/(fx - F(x0));
	    	x0 = x1;
			x1 = x1 - fx*s;
	    }
	    return x1;
    }
};
#endif
